!p.multi=[0,2,1]
!p.charsize=1.5

; Stratification code for ttauri

;-----------------------------
; Parameters
;-----------------------------
l = 188                    ;Radial steps
rr = 8.52e+08              ;Star`s radius in meters
b = 0.15                   ;Relative to the radius    
rb = b * rr                ;Bottom of domain 
t = 0.97                   ;Relative to the radius   
rt = t * rr                ;Top of domain  
irb = 1. / rb              ;Inverse of the bottom radius
Msol = 1.988e+30           ;Solar mass
dr = (rt - rb)/double(l)
r = rb + indgen(l)*dr      
M = 0.7 * Msol             ;Star`s mass
G = 6.674e-11              ;Gravitational constant
g = G*M / rr^2)  ;Gravity varying linearly with radius
r_ = indgen(l-1)*dr        ;Radius relative to the bottom of domain 
;------Bottom values:
gb = g[0]                  ;Gravity
rhob = 3.e+03              ;Density
Tb =  5.060e+06            ;Temperature
thetab =  5.060e+06        ;Potencial Temperature

Rg = 13732.0
cp = 2.5 * Rg
cap = Rg / cp
icap = 1. / cap

ss = gb / (cp * thetab)
rho = rhob * (1. - (ss * r_) / (1. + r_ * irb))^(icap - 1)

st = -3. / (thetab * (l - 1) * dr)
theta = thetab * (1. + st * dr)  

plot, r/rr, rho 
plot, r/rr, the, yr = [min(theta),max(theta)]
!p.multi=0
end
